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Abstract 


The  paper  describes  the  major  aspects  of  modeling  engineering  problems 
of  elastomechanics.  It  shows  various  aspects  and  results  on  a  set  of 
illustrative  examples  of  2  and  3  dimensional  problems. 
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1. 


Introduction 


The  aim  of  computational  analysis  is  to  describe  and  reliably  predict 
physical  phenomena  of  interest.  In  the  engineering  sciences  the  primary  aim 
is  usually  to  design  tools  which  operate  SAFELY  under  certain  (mechanical) 
conditions,  in  certain  environments  and  for  a  certain  period  of  time. 

By  computational  analysis,  ONLY  mathematical  problems  and  NOT  the  reality 
can  be  analyzed.  The  mathematical  problem  TRANSFORMS  given  input  data  into 
information  which  is  of  direct  interest  and  does  not  add  anything  new  (in 
fact  its  loses  some  information), 

The  aim  of  computation  is  to  reliably  obtain  certain  information  in  the 
range  of  an  admissible  tolerance  so  that  it  is  not  unduly  influenced  by  the 
computational  procedure  used. 

The  formulation  of  the  mathematical  problem  is  usually  the  most  crucial 
part  of  the  analysis.  Because  of  the  complexity  of  engineering  analysis  and 
uncertainties  in  the  available  information,  the  formulation  of  the  mathe¬ 
matical  problem  is  often  directly  or  indirectly  stipulated  in  the  design 
codes  and  often  (af  least  in  pants)  it  is  also  influenced  by  the  particulan 
(company)  engineering  practices.  These  codes  change  with  time  and  express  the 
experiences  with  the  technology  used.  As  a  typical  exajnple  we  mention  the 
design  code  (USAF-MIL-A-83444)  used  in  aircraft  components.  It  is  based  on 
the  principle  of  "non- inspect able  slow  crack  growth"  which  should  meet  the 
following  demands 

a)  the  life  of  the  component  should  exceed  two  design  times 

b)  the  residual  strength  of  the  component  should,  after  being  in  service 
two  design  life  time,  exceed  maximal  load  acting  on  the  component  by  a  factor 
of,  say,  9/8. 
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These  or  similar  principles  and  other  considerations  (for  example,  un¬ 
certainties  in  the  input  information)  lead  to  the  precise  formulation  of  the 
mathematical  problem  and  defining  the  data  which  have  to  be  obtained  as  well 
as  to  the  admissible  accuracy  with  which  they  have  to  be  determined. 

The  basic  flow  chart  of  an  engineering  computational  analysis  is  shown 
in  Fig.  1. 


1.  PHYSICAL  PROBLEM 
AND  CRITERIA 


2.  BASIC  MATHEMATICAL 
PROBLEM 


3.  SIMPLIFIED  MATHEMATICAL 
PROBLEM 

AND 

ANALYSIS  OF  THE  ERRORS 
CAUSED  BY  THE  SIMPLIFICATION 


4.  NUMERICAL  TREATMENT 
AND 

ANALYSIS  OF  THE  ERRORS 
CAUSED  BY  NUMERICAL 
TREATMENT 


5  PHYSICAL  CONCLUSIONS 
OR  ENGINEERING  DECISION 


Fig. 1.  The  flow  chart  of  computational  analysis. 
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Usually  in  practice  loops  are  present  in  the  flow  chart. 

The  "reality"  is  associated  with  (1).  The  engineering  analysis  of  the 
problem,  the  aims  of  the  analysis,  and  the  assessement  of  the  quality  of  the 
available  data.  etc.  then  yield  the  precisely  formulated  mathematical  problem 
(model)  (2).  This  model  is  to  be  understood  as  a  "higher"  model  which  is 
identified  with  reality.  Nevertheless  we  solve  usually  only  a  simplified 
problem  (3)  and  reliability  of  its  solution  is  judged  typically  in  comparison 
with  (2).  In  (4)  we  solve  numerically  problem  (3)  and  the  reliability  and 
the  error  of  the  numerical  solution  is  related  to  the  (exact)  solution  of  (3) 
(and  not  (2)  or  ( 1 ) ) . 

Let  us  underline  that  basic  and  simplified  mathematical  problem  has  to 
have  reasonable  mathematical  properties,  for  example  the  existence  of  a 
solution.  The  existence  of  the  solution  of  the  "real"  problem  does  not 
necessarily  mean  that  the  solution  of  the  mathematical  problem  exists  too. 

This  is  because  of  the  simplification  which  enters  into  the  formulation  of  the 
mathematical  problem.  Also,  if  the  numerical  algorithm  provides  numerical 
results  (possibly  reasonable  looking),  it  does  not  mean  that  the  solution  of 
the  mathematical  problem  necessarily  exists,  (because  convergence  has  not  to 
occur  etc. )  Obviously,  theoretical  analysis  of  the  mathematical  problem  and 
comparison  of  its  properties  with  the  expected  properties  of  the  reality  is 
essential  part  of  the  reliability  of  the  model. 

In  all  stages  we  have  to  relate  the  numerically  obtained  solution  to  the 
exact  solution  of  a  mathematical  problem.  The  agreement  with  reality,  for 
example  with  experiments,  is  then  related  solely  to  the  formulation  of  the 
basic  (or  possibly  simplified)  problem.  It  is  essential  that  the  errors  of 
numerical  solutions  are  completely  under  control  so  that  the  exact  solution  of 
the  mathematical  problem  is  essentially  achieved  and  a  possible  disagreement 
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with  experiments  is  related  ONLY  to  the  mathematical  problem  itself. 

In  this  paper  we  will  show  various  concrete  examples  to  illustrate  the 
basic  ideas  and  results.  All  the  computations  in  this  paper  has  been  made  by 
the  h-p  version  of  the  finite  element  method  by  the  code  PROBE  (McNeil 
Schwendler-Noetic)  and  STRIPE  (Aeronautical  Institute  of  Sweden).  These  codes 
have  various  error  checks  so  that  the  numerical  results  presented  here  can  be 
assumed  to  be  exact  in  the  range  of  accuracy  needed  for  the  model  conclusions. 


2)  Problem  of  the  cantilever  beam 


Consider  a  problem  of  a  simply  supported  cantilever  beam  shown  in 


Fig. 2. la 
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Fig. 2.1  The  simply  supported  cantilever  beam 

Let  the  basic  mathematical  problem  be  the  problem  of  two  dimensional  linear 

elasticity  (plane  strain)  for  isotropic  homogeneous  material.  The  basic 

unknowns,  the  displacement  u,v,  satisfy  the  usual  Lam6-Navier  equations  of 

elasticity.  As  the  boundary  conditions  shown  in  the  Fig. 2. la.  we  Impose  on 

BC  :  T  =  p,  T  =  0,  on  AB,  u  =  v  *  0,  on  CD,  DE,  FA  ;  T^  =  T  =  0  and  on 
y  x  3 

EFCwhere  &  «  d)  T  =  0,  v  =  b(x-L'  )  where  b  is  such  that  J*  T  (x-L'  )dx  =  0 
x  rr  y 
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By  respect  to  we  denoted  the  tract ior.r.  The  problem  is  a  model 
of  the  simply  supported  cantilever  beam.  The  (weak)  solution  which  has  finite 
energy  exists  and  is  uniquely  determined. 

Let  us  now  consider  the  simplified  problem  when  A  =  0  and  when  v  =  0 
at  the  point  G  (see  Fig.  lb)  is  prescribed  instead  of  the  more  complicated 
boundary  condition  of  the  basic  problem.  Then  it  Is  possible  to  show  that  the 
unique  (weak)  solution  of  the  simplified  problem  is  the  same  as  the  one  when 
the  condition  v  =  0  at  G  is  not  present.  We  see  that  the  solution  of  the 
simplified  problem  is  unacceptable.  The  reason  for  it  is  that  the  displace¬ 
ment  under  concentrated  load  is  infinite  (the  Bousinesque  solution).  Although 
the  solutions  of  the  basic  problem  converges  to  the  simplified  one  as  A— >0, 
the  convergence  is  very  slow  and  so  it  is  inadmissible  to  consider  this 
limiting  case  instead  of  the  original  one. 

We  remark  that  the  point  support  is  standardly  used  in  finite  element 
computations,  i.e.  the  simplified  problem  is  often  numerically  solved.  Hence 
the  error  of  the  finite  element  solution  is  very  large  because  for  a  mesh  not 
extremely  refined  the  cantilever  beam  solution  without  support,  i.e.  the  exact 
solution  is  not  obtained.  It  is  possible  to  show  that  the  finite  element 
solution  converges  to  the  solution  for  the  beam  without  support  as  the  mesh 
size  converges  to  zero.  Hence  the  solution  obtained  in  practice  is  mesh 
dependent.  This  is  of  course  completely  undesirable.  It  is  necessary  to 
mention  that  the  finite  element  meshes  used  in  practice  (if  not  adaptively 
constructed)  are  crude  and  the  FE  solution  does  not  show  the  mentioned  effect. 
For  some  numerical  analysis  and  computation  we  refer  to  [3). 
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3)  The  problem  of  the  built  in  plate  (beam). 


Let  us  consider  the  classical  problem  of  an  infinite  plate  problem  (in  2 
dimension)  which  can  be  formulated  as  two  dimensional  (plane  strain)  problem 
in  the  coordinates  xy.  The  scheme  of  a  concrete  example  is  shown  in  the 
Fig. 3.1  where  we  assume  built  in  (clamped)  boundary  conditions. 
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Fig. 3. 1  The  scheme  of  the  considered  plate. 


We  define  once  more  the  basic  problem  as  the  two  dimensional  elasticity 
problem  with  v  -  0.3  (where  we  denoted  by  v  the  "Poisson  ratio”)  and  with 

7 

the  modulus  of  elasticity  E  =  3. 10  ,  The  built-in  (clamped)  boundary 
condition  is  modeled  by  u  =  v  =  0, 

The  solution  exists,  and  is  unique.  Assume  that  the  aim  of  the  analysis 
are  the  stresses  in  the  areas  A  and  B  shown  in  Fig.  3.  1.  Let  us  further 
distinguish  2  cases  for  the  data  of  interest 

a)  the  bending  moment  and  the  shear  force 

b)  maximal  stresses  and  the  stress  distribution  through  a  cross  section. 


3.1.  The  problem  of  the  boundary  conditions. 

It  is  obvious  that  the  modeling  of  the  clamped  end  is  an  idealization. 
In  reality  the  support  is  obviously  more  complex.  Hence  we  have  an 
uncertainty  in  the  formulation  of  the  boundary  condition  and  the  problem 
shown  in  Fig. 3.1  can  be  understood  as  the  simplified  one.  To  analyze  this 
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problem  consider  a  few  configurations  which  could  be  expected  to  lead  the  same 
simplified  p -oblem.  They  are  shown  in  Fig. 3. 2. 


Fig. 3. 2  The  scheme  of  various  boundary  conditions. 

In  the  case  (e)  we  model  the  clamped  end  as  the  elastically  built-in  end. 

The  boundary  corditions  are  then 

T  =  -cv 
x 

T  =  -cw 

y 

where  c  =  10^. 

The  case  (f)  depicts  a  still  further  simplified  problem  based  on  the 
Klrchhoff  beam  theory  (strength  of  material  approach).  In  this  case  the 
value  of  the  bending  moment  in  the  center  is 


(3.  1 ) 
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in  the 


Assume  first  that  we  are  interested  in  the  maximal  sti ess  < r 

x 

center  of  the  beam  (area  A  in  Fig.  3.  1).  We  get  then  in  the  case  f  :  = 

100. 

We  let  the  stress  cr^  =  100.23  in  the  case  d)  (i.e.  the  case  shown 
in  Fig. 3.1)  be  the  "exact"  solution  to  which  we  will  compare  all  others. 
Table  3.  1  shows  the  results. 

Table  3.1.  The  stresses  at  the  center  for  various  models. 


Table  31  shows  that  simplification  of  either  problem  a,b,c  leads  to  the  error 
about  10%  while  the  simplification  of  the  case  d  by  Kirchhoff  hypotheses 
leads  to  error  of  0.2%. 

We  have  considered  the  special  case  for  the  ratio  d/L  =  If  d— >0 

(for  fixed  L),  the  relative  difference  between  the  cases  a,b,c,  and  d  goes 
to  zero.  We  have  then 
In  the  case  a): 


(3.2) 


aad(d)  ’ 


Cjd  +  higher  order  terms. 
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shows  that  in  fact  a^Cd)  =  with  high  accuracy. 


Table  3.2  .  The  relative  error  ^(d)  of  the  case  a  with  respect 
to  the  case  d 


d/L 

a  (d) 
ad 

1/20 

8 . 2% 

1/50 

3.  2% 

Let  us  now  consider  simplification  of  the  problem  d  to  the  problem  f. 
We  have  then 


*df(d) 


+  higher  order  terms. 


We  see  that  the  error  of  the  modeling  of  the  boundary  conditions 
is  much  more  significant  than  the  error  of  the  simplification  leading  to  the 
Kirchhoff  (strength  of  material)  solution. 

So  fair  we  have  computed  the  maximum  of  the  stress  in  the  area  A. 

Here  the  stress  is  very  accurately  linear  through  the  cross  section  and  hence 
when  the  interest  is  in  the  moments,  the  relative  errors  mentioned  in  the 
table  3.1  and  3.2  hold  too. 

We  reported  the  stress  in  the  area  A  and  have  seen  that  the  sensitivity 
to  the  boundary  conditions  is  of  order  10%. 

In  the  area  B  the  differences  in  the  stresses  are  much  larger.  The 
stresses  are  singular  (the  singular  behavior  of  the  solution  will  be  discussed 
in  the  next  section).  Here  we  report  in  the  Fig. 3. 3  a,b,c,  the  stress  in  the 
cross  section  in  the  distance  x  =  d/100  from  the  boundary.  We  clearly  see 
that  the  differences  between  the  mentioned  cases  are  significant.  On  the 
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Fig.  3.3  The  stresses  in  the  area  B. 


4)  The  singularity  problem  and  zooming  principles. 

4.  1  The  problem  in  2  dimensions. 

Let  us  consider  the  linear  elasticity  problem  on  a  polygon  domain  Q 

with  the  boundary  consisting  of  straight  segments  i  =  1 . n  and 

vertices  A^.i  =  lt...,n.  By  we  denote  the  internal  angles.  Let  us 
assume  that  we  are  dealing  with  a  homogeneous  isotropic  material  and  that  on 
every  segment  are  boundary  conditions  with  analytic  data  prescribed  and 
that  no  volume  forces  are  present.  Then  the  solution  is  analytic  on  nxUAj. 

The  solution  is  singular  in  the  neighborhood  of  the  vertices.  Let 
r^r^i  be  the  segments  meeting  in  the  vertex  A^.  Assume  for  simplicity 
that  the  boundary  condition  (of  standard  type)  are  homogeneous  on 
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Fig. 4.1  Scheme  of  the  angle 


Then  the  solution  of  the  problem  in  the  neighborhood  at  has  the  form 


(4.  1) 


u(x, y) 
v(x, y) 


n  A . 

- 1  v  1 


i=l 


<Pi  (0) 
(0) 


+  smoother  terms 


where  (r,0)  are  the  polar  coordinates  with  the  origin  in  as  shown  in 
Fig. 4.1.  are  real  a  complex  and  Re  Ai  +  1>ReAi  >0  and  <p ^  and  are 

smooth  functions  in  0.  If  A is  complex  then  we  use  real  and  imaginary 
parts  separately.  Coefficients  A^  and  <p^,  ^  are  given  by  the  geometry  (the 
angle  u^),  the  type  of  boundary  condition  on  +  1  and  the  material 

properties  (in  our  case  Poisson  ratio).  They  are  independent  of  the  solution. 
The  coefficients  Cj  depend  globally  on  the  solution  (except  for  special 
cases).  The  structure  of  the  solution  is  well  known  in  the  general  case,  for 
straight  and  curved  segments  T^,  and  for  general  (linear)  materials.  We  will 
not  go  into  details.  Here  we  refer  instead  to  [14],  [15],  [16],  [22],  In 
[21]  a  general  approach  (and  a  computer  code)  for  computation  of  A^,  <p  ,  \fi^ 
for  anisotropic  and  nonhomogeneous  materials  is  given.  Very  often  the 
coefficients  c^  called  stress  intensity  factors  (together  with  A^,  <p^,  \(i^,  ) 
are  the  main  aims  of  the  analysis  (in  2  and  3  dimensional  settings).  This  is, 
for  example,  in  the  case  of  design  based  on  the  earlier  mentioned  design  code 
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(USAF-MIL-A-83444) .  The  stresses  in  the  neighborhood  of  A^  are  unbounded 
when  ReA^  <  1  provided  that  c^*0. 

Because  in  finite  element  computations  the  stresses  are  always  finite, 
the  character  of  the  computed  stresses  can  be  misleading.  If  A^  <  1  then 
practically  always  (except  in  symmetric  cases)  c^*0  although  it  can  be 
relatively  small.  Then  the  large  stress  can  be  confined  to  small  area  only 
(we  will  see  an  illustration  in  a  three  dimensional  problem  in  the  next 
section. )  Reliable  computational  analysis  always  requires  to  compute  these 
stress  intensity  factors.  (For  methods  of  reliable  computation  of  the 
singular  behavior  around  the  corners  see  references  [1],  [25]). 

Let  us  now  relate  (4.1)  to  the  zooming  principle.  To  this  end  let 

_  x  y  r 

?  =  £•  v  =  K  and 

U(€,T?)  =  u(#c£,  K7j) 

V(C.TJ)  =  v(k£,  KTj ) 


be  the  zoomed  solution.  Then  we  obviously  have 


(4.2) 


U(?.T)) 


•  ^ 

•  l  C  F  U) 

i=l 


.  Ve) 


and  hence  functions  p^Vj(0)»  p^*^j(®)  are  the  parts  of  the  zoomed  solution 
with  the  zooming  parameter  k.  Here  it  was  characteristic  that  the  Infinite 
sector  was  invariant  with  respect  to  the  zooming. 

Let  us  consider  now  another  case,  namely  the  case  shown  in  Fig.  4.2.  We 
can  zoom  (Fig.  4.2b)  the  solution  at  A  (see  Fig.  4.2a)  and  get  the  solution 
in  the  form  (4.1)  (resp.  (4.2))  in  the  same  way  as  before.  The  other  possi¬ 
bility  is  to  zoom  the  solution  with  respect  to  the  parameter  d  as  shown  in 
Fig.  4.2b.  Then  up  to  rigid  body  motion  the  first  term  in  the  zoomed  solution 
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has  the  form 


(4.  3) 


u(c.-n) 

V(C.H) 


=  d-1M 


+  T  J 


*2(€.-o) 

l  02(€.t») 


Here  M  and  T  is  the  moment  and  shear  force  at  the  end  of  the  beam.  The 

function  «i>  i  =  1,2  are  the  functions  defined  on  the  domain  shown  in 
Fig.  4.2b. 


DETAIL  a 


'1 


DETAIL  b 

.L 


T 


Fig. 4. 2.  The  solution  on  the  zoomed  domain. 


Coefficient  M  and  T  are  analogous  to  the  stress  intensity  factors 
introduced  earlier.  The  zooming  principles  can  be  used  in  many  cases  when 
the  corner  singularity  interferes  as  in  the  case  shown  in  Fig.  4.2.  For  more 
details  we  refer  to  [61. 


4.3.  The  problem  in  3  dimensions 

In  3  dimensions  we  will  consider  a  polyhedron  instead  of  a  polygon  and 
the  problem  becomes  more  complex.  Once  more  we  will  consider  the  case  of 
isotropic  homogeneous  maeterial.  Along  the  edges  the  solution  is  singular  in 
the  direction  which  Is  perpendicular  to  the  edges  and  is  smooth  along  the 
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edges.  Assuming  that  the  (straight)  edge  is  along  the  axis  z,  the  singular 
terms  are  of  the  form 


u(x, y, z) 
v(x, y, z) 
w(x, y, z) 


(z)r*1 


<P\  (0) 

*h(e) 

<i(0> 


where  (r,0,z)  are  the  cylindrical  coordinates.  Functions  31,6 

smooth  in  0.  There  are,  in  contrast  to  the  two  dimensional  case,  two  kinds 
of  singular  function.  For  the  first  kind  coefficients  A^  and  functions  y  , 
i />.  are  as  in  2  dimensions  and  =  0  (they  are  sometimes  called  bending 
singularities).  In  the  second  kind  we  have  <p^  =  ^  =  0  (and  £^*0)  (they 
are  sometimes  called  torsion  singularities).  For  details  we  refer  to  [15], 
[22],  The  function  c^(z)  are  the  stress  intensity  functions  which  are  smooth 
in  z  (except  the  neighborhood  of  the  verticles). 

In  addition  to  the  edge  type  singularity,  we  have  a  vertex  singularity. 

Here  the  singular  terms  have  the  form 


u(x, y, z) 
v(x, y, z) 
w(x, y, z) 


(i) 


(Pi(e.H) 

^i(e.s) 

Ci(e.s) 


where  (R.0.H)  are  the  spherical  coordinates.  Functions  have 

singular  behavior  in  the  neighborhood  of  being  the  coordinates 

of  the  edges.  For  details  once  more  see  [15],  [22]. 

The  coefficient  A^,  and  functions  depend  on  the  geometry, 

boundary  conditions  and  material  properties  but  not  on  the  solution.  The 
coefficients  Cj  are  the  analog  to  the  stress  intensity  factors  and  depend 
(except  special  case)  globally  on  the  solution.  The  relation  between  A^ 
and  A^  governs  the  singular  behavior  of  the  edge  intensity  factor  function 
in  the  neighborhood  of  the  vertex.  The  unboundedness  of  the  stresses  in  the 
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neighborhood  of  the  vertex  take  place  when  A  <1  provided  the  coefficient 


0^*0.  Because  of  global  dependence  of  c^  on  the  solution,  c^*0 
practically  always  (except  possible  in  special  case  of  symmetries).  Neverthe¬ 
less  cam  be  small  and  large  stresses  cam  be  confined  to  a  small  area. 

Then  the  usual  finite  element  solution  could  indicate  completely  wrong 
behavior.  Hence  computation  of  c^  is  a  necessity  to  obtain  reliable  results 


Fig. 4. 3.  The  3  dimensional  domain 


TORSION  INTENSITY 
FACTOR  FUNCTION 


BENOING  INTENSITY 

factor  function 


Fig. 4. 4.  Location  of  the  rays 


Let  us  show  now  two  typical  examples.  Fig. 4. 3  shows  a  three  dimensional 
domain  with  imposed  boundary  conditions.  Along  the  marked  edge  shown  in  Fig. 
4.4  the  stress  intensity  factor  functions  are  present.  In  the  neighborhood 
of  A  solution  has  vertex  type  of  singularity.  In  Fig. 4. 4  we  show  the  rays 
where  the  stresses  are  depicted.  The  first  two  coefficients  A.  are  given  in 
Fig. 4. 5a  for  v  =  0.0  and  v  =  0.3.  The  stresses  on  the  rays  will  have  then 

the  form 

A(1)  1  A(2)-l 

a-  =  c  R  ♦  CLR  +  higher  order  terms. 

1  2 
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R.  and  relate  to  the  need  to  compute  two  stress  singular  functions. 

Fig .4.5  ab  depict  the  stresses  on  different  rays.  The  Fig.4.5ab  show  well 
also  the  scales  where  the  large  stresses  will  appear  (which  is  of  the  order  of 
1/100  of  the  thickness).  We  also  see  that  the  value  of  v  does  not  signifi- 
ntly  influence  the  behavior  of  the  solution.  Nevertheless  this  is  not  always 
the  case  as  we  will  see  in  the  next  example. 

The  second  problem  (which  was  suggested  by  K. J.  Bathe,  (see  also  [9])  is 
depicted  in  Fig.  4.6. 


CLAMPED 


Fig.  4.6  The  Bathe’s  problem 

The  boundary  conditions  are  shown  in  the  figure.  When  nothing  is  explicitly 
stated  then  the  associated  tractions  are  zero.  Let  us  consider  now  the 
behavior  of  the  solutions  at  the  edge  I-A:  for  v  =  0.0  and  i>  =  0.3. 

The  values  of  the  constants  are  the  same  as  in  the  previous  example. 

Fig. 4. 7  shows  now  the  stress  <r at  the  edge  I-A.  We  see  here  drastic 
difference  between  the  stress  behavior  for  v  =  0.0  and  v  =  0.3.  The 
standard  finite  element  analysis  will  lead  to  the  conclusion  that  a is 
bounded  on  the  edge  I-A  for  v  =  0.3  while  for  v  =  0.0  is  large.  This 
conclusion  is  of  course  completely  wrong.  Let  us  mention  that  the  strength  of 
material  solution  predicts  <r  =  -280  for  x  =  -d  and  z  =  0.  This  value 
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is  approximately  achieved  for  v  -  0.0  for  any  y  (the  problem  here  is 


•400  - 
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-300  | 
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- i/d - ► 

Fig. 4. 7  The  stress  < on  I-A  for  v  =  0.0  and  v  =  0.3 

* 

y-independent ) .  For  v  -  0.3,  it  is  approximately  achieved  for  y > 7d  . 

We  see  that  for  a  reliable  conclusion  about  the  stresses  in  the  neighborhood 
of  the  corners  and  edges  the  computation  of  the  stress  intensity  factor  is 
essent ial  and  any  program  should  always  have  to  be  able  to  provide  them 
(STRIPE  provides  them).  For  more  about  the  problem  of  reliable  computations 
in  solid  mechanics  and  detailed  analysis  of  Bathe  problem  and  engineering  we 
refer  to  [ 1  ] . 

5.  The  plate  problem 

The  plate  (and  shell)  problem  is  a  basic  problem  in  engineering.  As  the 
basic  mathematical  problem  we  will  understand  the  three  dimensional  problem 
of  elasticity  on  the  thin  domain 

Analysis  of  the  problems  of  this  section  has  been  made  by  programe  STRIPE  by 
Dr.B.  Andersson.  For  more  see  also  [ 1 ) . 
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-  (x,y,z|(x,y)  I  z  |  <  d/2} . 


We  will  assume  isotropic  homogeneous  material. 

5.1  The  problem  of  the  derivation  of  the  simplified  formulation. 

There  are  very  many  formulations  of  the  plate  problems.  For  the  survey 
see  [20]  [23],  The  major  simplified  formulations  are  the  Kirchhoff  [K]  and 
Reissner-Mindl in  (RM)  formulations.  Many  results  describe  the  asymptotic 
behavior  of  the  solution  of  3  dimensional  (basic)  formulation  as  d — >0. 

These  results  show  that  (for  example  in  the  energy  norm)  the  3D  and  Reissner- 
Mindlin  solutions  converge  to  the  Kirchhoff  solution  as  d — >0,  see  eg. 

[12],  [19].  For  a  detailed  study  of  the  asymptotic  behavior  of  Reissner- 
Mindl  in  problem  we  refer  to  [2]. 

Further  there  are  generalized  models  based  on  the  projection  in  the 
energy  on  the  space  of  functions  of  the  form  (Kantorovich  method) 

n  , 

(5.1)  U(x,y,z)  *  £  ¥>k(x,y)zK 

k=0 

n  , 

V(x.y.z)  *  £  t/»k(x, y)zK 

k=0 

m  ir 

W(x,y,z)  *  £  £k(x,y)z  . 

k=0 

In  the  case  v  =  0,  the  case  n  =  1  with  <Pq  -  <I>q  -  0  and  m  =  0  leads 
to  RM  model.  For  v  > 0,  n  =  1,  m  =  2,  s  8  s  0  leads  to  RM  model 

with  singular  perturbations.  Usually  m  =  n  ♦  1  is  taken,  which  guarantees 
the  proper  asymptotic  rate  of  convergence.  The  model  based  (5.1)  will  be 
called  n-m  model.  In  general  n  and  m  can  be  different  in  different  parts 
of  u.  For  n,m — »»  the  solution  of  (n-m)  model  converges  to  3  dimensional 
solution  (in  the  energy  norm). 
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The  form  (5. 1)  leads  to  a  hierarchic  family  of  models.  It  depends  on 
particular  choice  of  the  function  in  z.  In  (5.1)  polynomials  have  been 
used.  This  choice  is  optimal  in  an  asymptotic  way  when  d — >0  (see  eg. 
[12],  [26]).  Other  optimal  choices  can  be  considered  too  [24). 

The  error  of  the  various  models  have  to  be  judged  in  the  relation  to 
the  3D  solution  and  which  data  are  of  interest. 

5.2  The  problem  of  the  rhombic  simply  supported  plate. 

Let  us  consider  the  plate  shown  in  Fig. 5.1.  The  simple  support  car  be 
formulated  as 

a)  hard  simple  support 

b)  soft  simple  support 


Fig.  5. 1  The  rhombic  plate 

In  the  case  a)  we  assume  on  the  lateral  sides  w  =  0  and  u^  =  0  where  by 

u^  we  denote  the  displacement  in  the  directions  of  the  tangent  to  the  boundary 

of  w  (and  hence  T.  *0),  u  is  free  (i.e.  T  =0). 

t  n  n 

In  the  case  b)  the  only  constraint  on  the  lateral  side  is  w  -  0  (and  hence 

T^  =  T^  =  0).  The  K-model  cannot  distinguish  between  these  two  supports. 

Few  problems  now  appear 

i)  How  much  do  the  solutions  for  the  two  models  of  support  differ? 
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ii)  Which  support  does  the  K-model  describe? 
iii)  How  accurate  (with  respect  to  the  3D  model)  are  the  K  and  RM  models? 

Answers  to  these  questions  depend  strongly  on  how  we  measure  the  error.  It  is 
well  known  that  the  major  difference  between  K  and  RM  model  is  in  the 
boundary  layer.  For  a  detailed  analysis  in  the  case  when  the  boundary  of  w 
is  smooth  we  refer  to  [2] 

If  we  take  the  energy  norm  measure  of  the  error  and  the  load  uniform  on 
the  upper  side  of  the  rhombic  plate  the  relative  error  in  %  for  the  K-model, 
and  the  m  =  n  =  2  model  (which  leads  to  the  slightly  smaller  error  the  RM 
model)  and  for  the  soft  simple  support  is  given  in  Table  5.1  (3D  solution  is 
taken  as  exact ) . 


Table  5.  1  The  relative  energy  norm  error  of  K  and  m  =  m  =  2  model  for 
soft  simple  support  in  V. 


v  = 

0.0 

v  = 

0.  3 

d  = 

0.  1 

d  = 

0.01 

d  = 

0.  1 

d  = 

0.01 

K 

(2,2) 

K 

(2,2) 

K 

(2,2) 

K 

(2,2) 

CO 

o 

o 

39.56 

12.57 

11.87 

3.  50 

34.  52 

11.  18 

9.  88 

2.94 

o 

o 

00 

39.91 

12.  59 

12.23 

3. 57 

60° 

42.  24 

12.72 

15.  46 

4.  14 

40° 

45.  43 

13.60 

20.  50 

4.  24 

30° 

48.27 

IS.  41 

22.66 

4.  34 

44.68 

15.03 

18.91 

3,  68 

Table  5.2  shows  for  a  =  90  and  v  -  0  the  error  for  the  hand  support  (the 
3D  solution  with  hard  support  is  taken  as  exact).  The  error  of  the  3D  solu¬ 
tion  with  the  hard  support  with  respect  to  the  soft  support  of  3D  formulation 
is  34.7%  for  d  =  0.  1  and  1 1 .  IV.  for  d  =  0.01  For  more  see  [19] 
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Table  5.2  The  relative  enery  norm  error  for  the  K  and  f2.2)  model 


for  hard 

support 

in  ’/.. 

Model 

d  =  0.  1 

d  =  0.01 

K 

20.  31 

2.03 

(2.2) 

8.22 

0.68 

The  difference  between  models  eg.  K,  RM,  n-m  and  3D  model  is  largest  in 

the  boundary  layer.  It  can  be  in  fact  very  large. 

To  illustrate  this,  we  consider  the  square  plate  (a  =  90,  v  =  0.0).  Let 

and  Q ^  be  the  shear  forces  on  the  line  x  =  0.5,  0<y<0.5  computed 

from  the  3  dimensional  solution  for  d  =  0.01  and  soft  support.  They  are 

shown  in  Fig.5.2ab.  Realizing  that  the  K-model  leads  to  Q  =0,  we  see  that 

yz 

K-model  is  unreliable  for  these  data  of  interest. 
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Fig. 5. 2  The  shear  force 


Q  and  Q  at 
xz  yz 


x=0.5  and  0<y<0.05 


The  K-model  approximates  relatively  well  the  hard  support  but  not  the 
soft  one.  In  the  larger  distance  from  the  boundary,  the  K-model  is  usually 
usable  also  for  the  soft  support.  Although  the  K-model  approximates  the  hard 
support  the  usual  approach  is  to  modify  the  reaction  by  the  derivative  of 

the  twist  moment  to  get  the  soft  support  reaction.  This  leads  to  a  reasonably 
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good  approximation  of  Qxz-  In  Fig. 5. 3  Q  for  soft  and  hard  support  with 
twist  moment  adjustment  is  shown.  Nevertheless  no  adjustment  of  K-model  could 
give  reasonable  values  of  Q 


5.3  The  reaction  Q 

xz 

There  is  a  significant  theoretical  difference  between  hard  and  soft 
sample  support. 

In  [6]  we  have  analyzed  the  behavior  of  the  solutions  on  a  regular  n-eck 
polygons  inscribed  in  the  unit  circle.  We  have  shown  that  for  hard 

simple  support  the  solution  in  converge  to  a  solution  on  the  circle  S 
but  which  paradoxically  is  not  the  solution  on  S  for  the  hard  simple 
support.  This  happens  for  K,  RM  and  3D  model.  In  contrast,  in  the  case  of 

the  soft  support,  such  paradox  does  not  occur.  This  shows  that  the 
mathematical  problem  of  the  hard  support  has  a  property  which  we  would  not 
expect  in  "reality"  and  hence  am  increased  precaution  ha^  to  be  given  when 
hard  support  model  is  used. 

52  The  problem  of  the  solution  singularity  of  the  plate  models. 

Let  us  once  more  consider  the  square  plate  (a  =  90°),  d  =  0.01,  v  =  0.3 
such  that  the  sides  0.5<x<0.5,  y  =  ±0.5  are  clamped  (u,v,w  =  0)  atnd 
other  two  sides  are  free.  Then  in  the  neighborhood  of  the  vertices  the 
solution  is  singular.  In  (111,  [27]  the  singular  behavior  of  the  RM  solution 
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is  analyzed.  This  analysis  shows  that  RM  solution  has  two  different 
characters  of  the  singularities.  Denote  by  (r,0)  the  polar  coordinates  with 
center  in  A  =  (0. 5,0.5).  Then  for  r  «  d  the  solution  for  example  the 
stresses  or  moments  have  the  form 

Mx  *  %  -  CrX'%M(e) 

while  for  r  ~  d,  the  singularity  is  as  for  the  Kirchhoff  model 

cr  s  CrAV^0) 

X 

Between  these  two  areas  there  is  a  transition  domain.  The  exponents  A 

KM 

\ satisfy  some  transcendental  equations.  In  our  case 

arm  =  -  °'241 

XK  =  +  0.0686  +  iO. 438 

(i.e.  the  stress  of  K-solution  is  oscillating).  In  Fig.  5.4  we  show  the 
stress  <r^  at  the  diagonal  of  the  plate  in  log  -  log  scale.  We  see  clearly 
that  both  types  of  singularities  occur.  Other  stresses  show  similar  behavior. 


Fig.  5.4  Stress  <r  of  the  RM  model. 

Finally  we  can  compute  the  character  of  the  moments  computed  from  the  3 
dimensional  solution.  Here  we  can  show  that 
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and  in  our  case 


A 


3D 


-0.289 


We  see  that  the  corner  behavior  is  different  for  these  3  models.  This 
difference  strongly  depends  on  the  geometry  of  the  plate  and  boundary 
conditions.  For  more  we  refer  to  [8]. 


The  problem  of  nonlinear  elasticity 

The  nonlinear  formulation  in  the  theory  of  elasticity  stems  from 
a)  nonlinear  geometry  as  large  displacements,  stresses,  etc. 


b)  nonlinear  constitutive  law. 


Here  we  will  address  some  questions  related  to  the  elasticity  assuming 
static  behavior  where  effects  of  velocity  etc.  can  be  neglected.  In  one 
dimensional  case,  given  the  strain  e(t)  -» < t < »  the  constitutive  law 
leads  to  the  stress  response  <r(t) 


(6.  1) 


<r  =  Ac 


where  A  is  an  operator  mapping  the  strains  into  the  space  of  stresses. 

In  the  3  dimensional  case,  c  and  cr  are  the  strain  and  stress  tensors, 
respectively. 

Usually  the  constitutive  law  in  2  and  3  dimensions  is  derived  from  the 
one  dimensional  law  by  applying  various  principles  as  Mises,  Tresca, 
Huber-Hencky,  etc.  In  one  dimension  mauiy  laws  were  proposed,  see  e.g.  [281. 
The  basic  laws  are  kinematic,  isotropic  hardening  and  others.  Recently  the 
formulation  by  Chaboche  [13]  has  become  popular. 

Mathematically  it  is  important  that  the  constitutive  law  is  such  that  it 
satisfies  conditions  which  guarantee  the  desirable  properties  of  the  mathe- 
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matical  problem  of  elasticity  where  it  is  used.  It  is  of  course  also 
important  that  there  is  not  a  large  difference  between  observed  and  predicted 
response  (based  on  the  constitutive  law  used). 

6.  1  Experimental  results 

Results  of  an  extensive  one  dimensional,  experimental  analysis  with  the 
aluminum  alloy  5454  in  the  H32  condition  are  reported  in  [18].  This  alloy  is 
produced  (under  the  same  commercial  mark)  by  different  manufacturers  and  is 
widely  used  in  engineering. 

The  analysis  in  [18]  is  based  on  the  fact  that  in  engineering  the 
material  is  taken  from  the  warehouse  and  at  best  the  experiments  for  selection 
of  the  proper  constitutive  law  can  be  made  on  samples  only  (statistical 
approach).  Hence  84  samples  have  been  taken  and  analyzed.  Among  others,  the 
main  questions  were  related  to 

a)  reproducibility  of  the  response 

b)  selection  of  the  constitutive  law. 

Two  classes  of  the  strain  were  considered 

i)  the  cyclic  periodic  strain  (which  is  usually  used  in  material 
science) . 

ii)  random  strain  which  is  more  realistic  in  applications. 

The  main  results  can  be  broadly  characterized  as  follows 

a)  The  reproducibility  factor  Q  for  the  random  strain  is  of  order 
-10-15%  where 


max | A( t )  -  B( t ) | 


Here  A(t)  and  B(t)  is  the  stress  response  of  two  different  sample  to  the 
same  random  strain.  For  cyclic  load  the  factor  Qc  is  of  order  7-10%. 
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b)  For  every  sample  and  particular  constitutive  law  mostly  used  in 
practice  (as  Chaboche,  kinematic,  Mroz,  etc.)  the  constants  for  the  best  fit 
were  computed.  Then  the  average  value  of  these  (84  samples)  were  computed  and 
using  these  constants  the  constitutive  law  the  factors  CR  (resp.  C  )  were 
analogous  to  QR(resp  Q  )  i.e.  we  define 


max  |  A  ( t )  -  B  ( t )  | 
_t _ 

max|A(t)  +  B(t) 


Here  A(t)  is  the  response  for  a  sample  and  B(t)  is  the  predicted  response 
based  on  the  average  constants.  For  the  best  law  (one  of  them  is  Chaboche  we 
get  CD  «  22-25*/.  and  C  *  16-18*/.. 

K  C 

For  the  best  fit  of  one  sample  we  get  C_  »  8%. 

R 

If  the  set  of  averaging  is  small  we  can  get  CD > 30%. 

R 

For  some  laws  (In  standard  use  in  FE  codes),  Co>  40-50*/.. 

R 


In  the  Table  6.1  we  show  the  il*llL  IMI^  norm  (in  psi)  and  relative  error. 
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l 

Kinematic 

1 

11850 

1 

3475 

32.8% 

17.6% 

(fd,  fq  is  the  label  of  the  sample,  Chaboche  and  Kinematic  means  response 
obtained  by  the  Chaboche  resp.  kinematic  law).  We  mention  that  for 


computational  purpose,  the  norm  11*11^  essential  (and  not 

00 


). 


As  an  illustration  we  show  in  Fig,  6.1,  6.2,  6.3,  the  value  for  A-B  (two 
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Fig. 6. 3  The  difference  between  sample  and  kinematic  law. 

[18]  analyzes  only  the  one  dimensional  problem.  It  is  possible  to  expect 
that  in  2  resp.  3  dimensional  setting  the  factors  will  be  larger. 

The  analysis  made  in  [18]  indicates 

a)  It  is  necessary  to  analyze  the  reliability  of  constitutive  laws 
derived  from  statistical  sampling. 

b)  It  is  necessary  to  analyze  random  and  not  cyclic  strains. 

c)  It  is  highly  desirable  to  develop  a  mathematical  theory  for  determin¬ 
ing  the  constitutive  law  based  on  the  (infinitely  dimensional)  identification 
problems  and  to  develop  a  strategies  for  optimal  selections  of  strains  for 
experiments.  This  is  especially  important  for  2  and  3  dimensional  settings. 

d)  It  seems  that  the  usual  elasticity  formulation  and  computation  based 
on  the  "average"  constitutive  law  cannot  give  reliable  results  and  other 
approaches  such  as  bracketting  have  to  be  developed,  see  also  here  [10]. 


6.2  Mathematical  formulation  of  Chaboche  law. 

We  have  seen  in  the  Section  6. 1  that  the  Chaboche  law  is  one  of  the 
laws  which  fits  best  the  data  for  the  single  sample.  Therefore  we  will 
discuss  it  here  in  more  detail. 

Although  in  [13]  the  law  is  formulated  in  an  incremental  way  related  to 

mechanical  interpretation,  it  can  be  cast  into  a  system  of  ODE  for  the  stress 

and  two  (internal)  parameter  functions,  ( <r( t ) ,  *(t),  R( t ) )  for  given  strain 

•  dc 

c(t).  In  what  follows  we  denote  e(t)  =  -rr-  etc.  The  Chaboche  law  is 

dt 

characterized  by  6  constants. 

We  have 

or  =  Ee,  <r(Q)  =  *(0)  =  R(0)  =  0 

X  =  0,  G^O)  =  Cy>  =  ~Cy 


(6.2) 


R  =  0 


eh  =  ° 
cl  =  0 


for  all  t € g,  where 


§  =  <t|e,(t)  <  e(t)  <  e.  (t),  or  e(t)  =  e.  (t) 

t  n  n 


and  c<0  or  e^(t)  =  e(t)  and  c  2:0) 


_  E[c(a-x(t))  +  b(Q-R(t))] 

^  =  cTa-*U) )  +  bT0~R( 1 5T"+  E  e(t) 


(6.3) 


_  E[c(a-x(  t ) )  Ifn 

Z(t)  '  cfi-iCt))'  ♦  bfQ-R(t)T+  E  G(t) 


ftm  _  Eb(Q-R(t) )  . 

R(t)  ‘  c(a-Z(t))  +  b(Q-R(t))  +  E  G(t) 


Gh  *  G 


G£  *  e  -  2E 
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for  all  t  €  P+  where 


P+  =  { 1 1  e  >  0  and  e 


Ch} 


and 


E[c(a-x(t) )  +  b(Q-R( t ) ) 
c(a+*(t) )  +  b(Q-R( t ) )  +  E 


R(t) 


tb(Q-R(t ) ) 

cta+z(t))  +  b(Q-R( t ) )  +  E 


e 


l 


for  all  t  €  P  ,  where 

P  =  { 1 1  e  <  0  and  e  =  e^> 

Chaboche  model  is  characterized  by  6  constants  with  a  physical  interpretation. 

a:  Kinematic  coefficient 
c:  Kinematic  exponent 
Q:  isotropic  exponent 
b:  isotropic  exponent 
g:  yield  strain 
ET:  elastic  modulus 


The  Chaboche  law  as  formulated  can  be  generalized  into  2D  and  3D  formu¬ 
lations.  Nevertheless  when  this  is  used  in  the  nonlinear  elasticity  equation 
problem,  some  desirable  mathematical  properties  of  the  problems  (as  for 
example,  the  existence  of  the  solution)  are  not  quaranteed.  Hence  another 
formulation  which  approximate  well  the  Chaboche  law  and  leads  to  the  desirable 
properties  of  the  mathematical  formulation  should  be  used. 
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6.3  A  proper  mathematical  formulation  of  the  constitutive  law. 

Let  us  outline  now  here  principles  and  family  of  constitutive  laws 
(called  gauge  method)  which  guarantee  good  properties  of  the  mathematical 
problems  based  on  them.  There  are  two  basic  (sufficient)  conditions  for  it. 

a)  Existence  and  convexity  of  the  yield  surface 

b)  the  normality  condition 

(These  conditions  are  related  to  the  Druckers  postulates). 

Let  ae!Rm  be  the  set  of  internal  parameters,  acAcRm,  A  being  a 
m  3 

convex  set  in  !R  .  Set  further  creR  (for  a  two  dimensional  problem). 

Now  we  will  formulate  the  law  with  the  help  of  the  yield  function.  To 

3 

this  end  let  F(<r,  a),  F  :  (R  xA— »R  be  given  so  that 

(6.8a)  F  is  convex  and 
(6.8b)  F(0, 0)  =  0 

(6.8c)  There  exist  constants  y,  F  such  that  0 < y <  |5  F||5  F|  < T 

a  a 

uniformly  on  the  set  { (<r,  a)  |  F(<r,  a)  =  z^>  for  some  z^. 

Then 


(6.9) 


(6. 10) 


where 


<r  =  Dc  if  teg 
a  =  0 

D6  F(5  F)TD 

_ or  or _ 

(5  F)T(5  F)  +  (5  F)T(5  F) 
a  a  <r  <r 

a  =  -  ( (6  F)T  -  (5  F)  )_1  ( (5  F)T<r)5  F  if  t€? 
a  a  <r  a 

g  =  { 1 1  F(<r,  a)  <  Zq  or 

F(a,a)  =  z  and  (5  FJ^irSO) 

0  a 
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7  =  { 1 1 F  ( <r ,  a )  =  zn  and  (5  F )  T<r  >  0 )  . 

U  <T 

The  Chaboche  model  could  be  cast  approximately  in  the  above  frame  using 

F(<r,a,6)  =  [maxlFj^  (<r,  a.  £) ,  F^o-.a.f?)] 

where 

(6.  10a)  Fj  (<r,  a.  j3)  =  (a-a^  +  a^fa-o^ )  +  a^O-gj  +  a^g-^  )  +  (<r-cr1 )  + 

(6.10b)  F2(<r,a.0)  =  bj (a-o^)2  +  a-ctg)  +  b^-^  )2  +  -  (<r-<r2 )  +  <2 

* 

and  [  ]  the  smoothing  the  operator  in  the  neighborhood  of  the  manifold 

F1l<r,a,0)  =  F2(<r,a,fi). 

Fig.  6.4a  shows  the  relation  between  e  and  <r  for  cyclic  (sinusoidal) 
strain  with  50  reversals  (25  periods)  for  the  constant  computed  out  of  the 
experimental  data  (averages  of  Chaboche  constants). 


Fig. 6. 4  The  relation  between  the  strain  e  and  the  stress  for  the  law  based 
on  (6. 10a, b) . 
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Fig.  6.4b  shows  the  results  where,  for  the  constant,  we  have  taken  the  mean 
minus  standard  deviation  (let  us  mention  that  the  correlation  for  these 
coefficients  (see  [18])  are  of  order  0.2).  Fig.  6.4c  shows  the  results  when 
computed  from  the  Chaboche  law  (see  section  6.2)  and  Fig.6.4d  shows  the 
experimental  results  for  one  sample. 

We  see  that  the  results  from  the  original  Chaboche  law  are  well 
approximated.  In  all  these  data  we  have  assumed  that  the  initial  data  which 
depend  on  past  history  have  been  known.  In  reality  they  are  not  known  which 
further  increases  the  uncertainties  in  the  available  information.  The 
derivation  of  2  and  3  dimensional  constitutive  law  leads  to  still  more 
uncertainties  because  not  enough  experiments  could  be  made.  In  [10]  the  2 
dimensional  constitutive  law  was  derived  as  the  limit  of  the  frame  made  out  of 
the  bars,  analogous  Cauchy’s  derivation  of  linear  elasticity. 

We  have  seen  that  the  computation  of  the  problem  of  elasticity  on  the 
assumption  of  the  knowledge  of  constitutive  law  without  respecting  the 
uncertainties  leads  to  unreliable  results. 

7.  A  posteriori  error  analysis  of  the  model. 

It  is  essential  to  make  a  posteriori  analysis  of  the  error  of  the 
solution  of  the  simpliefled  problem  when  only  the  data  from  this  simplified 
model  are  used.  This  cam  often  be  made  by  two  sided  energy  estimates.  For 
details  see  e.g.  [17]. 

7.1,  Estimate  of  the  error  of  geometry  idealization. 

Consider  the  problem  on  shown  in  Fig.  7.  1, 
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Fig,  7.1.  Scheme  of  the  problem  with  perturbed  boundary. 


Let  the  basic  mathematical  problem  be  the  1  inear  elasticity  on  0^  and 

the  simplified  problem  is  the  problem  on  £1^  using  E  =  1,  i>  =  0.3  (and 

r  =  0).  Then  by  finite  element  solution  we  get  two  stress  intensity  factors 

2  2 

(see  Section  4)  Cj  =  0.2157  10  ,  c^  =  0.5929  10  .  Then  the  estimate  using 

Cj,  c2  (and  the  form  of  the  singular  functions)  allows  to  compute  the  upper 

estimate  in  the  energy  norm  on  nQ  of  the  difference  between  the  solutions  on 

and  (see  [17]).  Table  7.1  gives  the  results  together  with  the  true 

error  obtained  by  tho  solution  on  0  . 

r 

Table  7.  1  The  estimate  and  true  error  of  the  geometry  ideal 


r 

Estimate 

True  error 

0.  1 

19.0% 

13.  2% 

0.01 

3.  5% 

3.0% 

We  see  good  effectiveness  of  the  estimate. 

7.2  Estimates  of  the  linearization 

Let  us  consider  once  more  the  problem  shown  on  Fig. 7.1  with  r  =  0.  If 
the  simplified  problem  will  be  understood  as  the  linear  problem  the  strains 
and  stresses  are  infinite  (see  Section  4).  If  we  would  consider  as  the  basic 
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problem  the  nonlinear  problem,  the  nonlinearity  will  occur  in  the  neighborhood 
of  the  corner.  We  will  assume  the  Hencky  model  of  the  nonlinear  elasticity 
with  the  governing  function 


*<o  =  \ 


<  for  0<<<C0 


with  s  € 


f0.  5-ff 
[  1-0  * 


1). 


Functions  is  the  function  describing  the  nonlinearity  of  the  material. 

g 

Let  us  consider  the  problem  depicted  in  Fig. 7.1  with  E  =  10  ,  v  =  0.3. 
Table  7.2  shows  the  upper  estimate  of  the  relative  error  mentioned  in  the 
energy  norm  for  various  values  of  For  more  details  we  refer  to  [  17 ] . 


J~t 

o 

II 

o 

o 

<0  =  0.001 

Estimate  of 

Estimate  of 

0 

7 

relative 

0 

7 

re  1  at  i  ve 

error 

error 

0.  1 

0.501 

2.76  10~4% 

0.9 

0.8 

2.05  10"5% 

0.  1 

0.5001 

2.79  10"4*/. 

0.9 

0.5 

4.06  10'5% 

0.  1 

0.50001 

2.79  10”4</. 

0.9 

0.02 

5.36  10~5% 

0.01 

0.501 

9.56  10”4'/. 

0.  5 

0.8 

1.  19  IO-4’/. 

0.  01 

0.5001 

9.62  10-4’/. 

0.5 

0.5 

2.55  10~4*/. 

0.01 

0.50001 

9.63  10~\ 

0,  5 

0.2 

3.89  10_4% 

0.001 

0.501 

3.03  io-3y. 

0.  3 

0.3 

1.84  10-4% 

0.001 

0. 5001 

3.07  io ~3y. 

0.  3 

0.5 

4.32  10-4’/. 

0.001 

0.50001 

3.08  10"3*/. 
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8. 


Cone  1  us  ions 


We  have  discussed  various  aspects  of  reliability  without  trying  to  define 
more  precisely  what  does  it  mean.  We  have  seen  that  tne  aim  is  to  get  desired 
data  with  an  assessment  of  their  accuracy.  More  precisely  we  are  aiming  to 
get  the  quantitative  bracketts  in  which  the  "true"  results  are.  These 
bracketts  then  express  the  uncertainty  of  the  results  caused  by  uncertainties 
in  the  input  data,  the  (simplified)  formulation,  the  discretization  etc. 

We  can  now  roughly  define  the  reliability  of  engineering  computations  in 
the  relation  to  reality. 

The  computational  results  furnished  with  the  bracketts  are  physically 
reliable  if  the  physically  observed  results  are  in  the  provided  bracketts. 

Analogously,  (more  precisely)  we  can  define  the  reliability  of  the  compu¬ 
tational  results  in  the  relation  of  mathematical  analysis. 

The  computational  results  furnished  with  the  bracketts  are  mathematically 
reliable  if  the  exact  data  of  the  basic  mathematical  problems  are  in  the 
provided  bracketts. 

We  have  seen  that  the  reliability  is  related  to  the  data  of  interest  and 
the  definitions  what  is  meant  by  accuracy  (e.g.  particular  norms)  etc. 

We  have  also  seen  that  the  mathematical  formulation  has  to  be  closely 
related  to  the  engineering  analysis  and  experimentation.  Without  it,  the 
"physical"  reliability  is  impossible  to  expect. 

The  mathematical  reliability,  i.e.  the  comparison  of  the  obtained  results 
with  the  exact  data  stemming  from  the  basic  mathematical  problem  is  always  (at 
least  in  principle)  possible.  This  comparison  and  its  bracketting  is  then  the 
main  goal  of  the  (mathematical)  computational  analysis. 

The  reliability  of  the  computational  analyses  has  many  features  and 
bring  out  many  unsolved  problems.  Nevertheless  there  are  already  today  many 
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ways  to  get  at  least  partial  quantitative  insight  into  reliability  of  computed 
results . 
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The  Laboratory  for  Numerical  analysis  is  an  integral  part  of  the 
Institute  for  Physical  Science  and  Technology  of  the  University  of  Maryland, 
under  the  general  administration  of  the  Director,  Institute  for  Physical 
Science  and  Technology.  It  has  the  following  goals: 

o  To  conduct  research  in  the  mathematical  theory  and  computational 

implementation  of  numerical  analysis  and  related  topics,  with  emphasis 
on  the  numerical  treatment  of  linear  and  nonlinear  differential  equa¬ 
tions  and  problems  in  linear  and  nonlinear  algebra. 

o  To  help  bridge  gaps  between  computational  directions  in  engineering, 
physics,  etc.,  and  those  in  the  mathematical  community. 

o  To  provide  a  limited  consulting  service  in  all  areas  of  numerical 
mathematics  to  the  University  as  a  whole,  and  also  to  government 
agencies  and  industries  in  the  State  of  Maryland  and  the  Washington 
Metropolitan  area. 

o  To  assist  with  the  education  of  numerical  analysts,  especially  at  the 
postdoctoral  level,  in  conjunction  with  the  Interdisciplinary  Applied 
Mathematics  Program  and  the  programs  of  the  Mathematics  and  Computer 
Science  Departments.  This  includes  active  collaboration  with  govern¬ 
ment  agencies  such  as  the  National  Bureau  of  Standards. 

o  To  be  an  international  center  of  study  and  research  for  foreign 

students  in  numerical  mathematics  who  are  supported  by  foreign  govern¬ 
ments  or  exchange  agencies  (Fulbright,  etc.) 


Further  information  may  be  obtained  from  Professor  I.  Babuska,  Chairman, 
Laboratory  for  Numerical  Analysis,  Institute  for  Physical  Science  and 
Technology,  University  of  Maryland,  College  Park,  Maryland  20742. 


